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Merging galaxy clusters such as the Bullet Cluster provide a powerful testing ground for indirect 
detection of dark matter. The spatial distribution of the dark matter is both directly measurable 
through gravitational leasing and substantially different from the distribution of potential astro- 
physical backgrounds. We propose to use this spatial information to identify the origin of indirect 
detection signals, and we show that even statistical excesses of a few sigma can be robustly tested 
for consistency—or inconsistency—with a dark matter source. For example, our methods, combined 
with already-existing observations of the Coma Cluster, would allow the 3.55 keV line to be tested 
for compatibility with a dark matter origin. We also discuss the optimal spatial reweighting of 
photons for indirect detection searches. The current discovery rate of merging galaxy clusters and 
associated leasing maps strongly motivates deep exposures in these dark matter targets for both 
current and upcoming indirect detection experiments in the X-ray and gamma-ray bands. 
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I. INTRODUCTION 

Over the course of the last century, astrophysical and 
cosmological observations have provided conclusive evi¬ 
dence of a non-luminous form of matter—dark matter 
(DM)—which comprises most of the matter content in 
our Universe. Signatures of the gravitational influence 
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of the DM energy density have been seen in the ro¬ 
tation curves of galaxies, orbital velocities of galaxies 
in galaxy clusters, strong and weak gravitational lens- 
ing, the cosmic microwave background, and large scale 
structure [T]. Nevertheless, the mass and possible non- 
gravitational couplings of the DM particle(s) remain un¬ 
known. There exist relatively weak model-independent 
upper and lower bounds on the DM mass, as well as 
a variety of constraints on its non-gravitational interac¬ 
tions from astrophysical observations, collider searches, 
direct detection experiments, and indirect detection ex¬ 
periments (see Refs. Hi] for reviews). 

The latter category of experiments searches for pho- 
toni[^ produced via decays or (semi-)annihilations of DM 
particles in astrophysical objects such as the galactic cen¬ 
ter (GC), dwarf spheroidal galaxies, or clusters of galax¬ 
ies. For example, a DM particle decaying to two photons 
would produce a monochromatic excess of photons with 
energy equal to half the DM mass—a line in the photon 
energy spectrum. A variety of experiments have been 
performed or are currently underway, particularly in the 
X-ray and gamma-ray energy bands. 

Most of the above searches have so far only yielded 
constraints, though several tantalizing anomalies have 
occasionally appeared in the data. Evidence for an in¬ 
tense line at 511 keV from the GC has accumulated over 
the last 30 years (see Ref. H for a review). Several 
groups have identified a broad excess of 1-10 GeV pho¬ 
tons originating in the GC in the Fermi-LAT data dDl- 
[T6| . Most recently, a line-shaped excess of 3.55 keV pho¬ 
tons has been identified in a spectrum of stacked galaxy 
clusters m in both XMM-Newton and Chandra data, 


^ High-energy neutrinos 0, or electron-positron pairs [8] are also 
interesting, but the methods in our paper are not relevant for 
these detection channels because of their lack of angular resolu¬ 
tion. 
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and independently in the Andromeda galaxy and the 
Perseus cluster [18] and the galactic center [19] with 
XMM-Newton data. 

A monochromatic photon line is usually considered to 
be a “smoking gun” signature of dark matter, but the 
energy spectrum of the signal may not be sufficient to 
confirm a dark matter origin of an excess. The positrons 
needed to fuel the bright 511 keV line in the GC may 
be provided by radioactive decays of unstable nuclei pro¬ 
duced in supernovae or massive stars, though a signifi¬ 
cant contribution from DM sources is not excluded [S]. 
While consistent with DM annihilation, the broad excess 
at 1-10 GeV in the galactic center could be due to other 
astrophysical sources such as unresolved millisecond pul¬ 
sars m- Lastly, the 3.55 keV line is consistent with, e.g., 
a 7.1 keV sterile neutrino decay, but could also arise from 
known emission lines in the same energy range |20] . 

When does a positive signal constitute a discovery? 
Additional input beyond the energy spectrum alone may 
be necessary to distinguish a DM signal from other astro- 
physical processes. The spatial distribution of the signal 
in the sky is the most obvious handle for this purpose. 
One can mask known point sources from the analysis 
(e.g. Sgr A* in the GC), or test whether an excess is 
consistent with DM annihilation or decay given a cer¬ 
tain model for the DM mass density p{r) in the Milky 
Way (MW). However, significant uncertainties arise in 
any quantitative analysis because the radial DM density 
profile in the inner regions of the MW has not been ex¬ 
perimentally measured and can only be estimated from 
Wbody simulations m- In addition, the spatial dis¬ 
tribution of possible background sources such as pulsars 
and supernova remnants in the MW remains poorly con¬ 
strained. Hence, the spatial distribution of an excess 
(as compared to the distribution of DM and background 
sources) in the MW cannot be considered a definitive 
test of a DM signal. The same is true for sources such 
as globular clusters and nearby dwarf spheroidal galax¬ 
ies (dSphs), whose DM content has only been measured 
indirectly (e.g. by the motion of their stars): the spa¬ 
tial profile of DM in these objects is typically an in¬ 
put for the calculation of their total mass. In galaxy 
clusters, however, the mass distribution of DM can be 
measured with weak gravitational lensing. Moreover, the 
major source of backgrounds in clusters—the intracluster 
medium (ICM)—emits X-ray radiation and therefore has 
a measurable spatial distribution. 

We advocate performing indirect detection 
studies in merging galaxy clusters, where the 
dark matter is physically separated from the 
baryons—and thus most backgrounds. The spatial 
separation between the DM and the bulk of the baryons 
can be exploited to construct powerful, quantitative tests 
of a potential DM signal. Although the GG and dSphs 
are expected to be much brighter sources of DM, and 
therefore more promising for the initial detection of a 
signal, merging clusters offer the unique possibility of 
comparing the spatial distribution of a signal with the 
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FIG. 1. (Color online) Weak gravitational lensing con¬ 
tours and X-ray luminosity maps of (a) the Bullet Clus¬ 
ter |22[ 1281 [2^ and (b) the Coma Cluster |26[I30 |. displaying 
the partial separation of the dark matter and the intraclus¬ 
ter medium. The lensing maps of (a) & (b) were smoothed 
with Gaussian kernels of cr = 20” & 1.7', respectively; con¬ 
tours signify integer multiples of signal-to-noise ratios (with 
S/N = 1 at the outermost contours). The X-ray density 
maps are proportional to the intensity of X-ray photons in 
the (a) 0.2-15 keV and (b) 0.2-2 keV energy band. Gridlines 
represent the binning (30” & 4', respectively) to be used in 
quantitative results in Sec. |IV[ 


known, distinctive distribution of DM. 

The first, spectacular evidence of a merger event caus¬ 
ing separation of dark and baryonic matter was observed 
in the galaxy cluster IE 0657-558, the Bullet Cluster [22] . 
Since then, lensing map reconstructions have identified 
many similar clusters. Major merger events typically 
lead to increased radio emission on the outskirts of galax¬ 
ies [23l |24] . Giant radio arcs are thus an easy “tag” of 
merging clusters, leading to the rapidly increasing avail¬ 
ability of weak lensing maps for this class of clusters. 
Most notably, the nearby A1656—the Coma Cluster— 
was discovered to have undergone a merger leading to a 
separation of the DM and ICM PSIBT] . In Table [ij we 
list relevant properties of the Bullet and Coma Clusters 
as well as several other merging clusters. We will use the 
Bullet and Coma Clusters as the two benchmark clusters 
of our paper; their lensing contours and X-ray luminosity 
maps are shown in Fig. 
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Z 

A0DM 

M200 

r200 


[Gp^ 


[arcmin] 

[h-Ho^^MQ] 

[h~^ Mpc] 

A1656 “Coma” [261I3T] 

0.10 

0.0231 

30 

1.9 

2.0 

IE 0657-558 “Bullet” [32] 

1.5 

0.296 

2.7 

2.‘E 

1.(0 

A3376 03 El 

0.20 

0.046 

icQ 

0.3 

1.1 

A520 “Train Wreck” 03 06] 

0.97 

0.199 

3 

0.7 

1.3 

A2163 07] 

0.99 

0.201 

1.6 

2.0 

1.8 

A1758 03 

1.4 

0.279 

1 

1.^ 

1.0 

A2744 “Pandora” [39ll40] 

1.6 

0.308 

2 

1.5 

1.7 

DLSCL J0916.2-k2951 “Musket Ball” 01] 

3.1 

0.533 

2.8 

2.2 

0.7 

MACS J0025.4-1222 02] 

3.4 

0.586 

1.3 

- 

- 


^ assuming Hq = 70kms ^ Mpc ^ {h = 0.7), Qa = 0.7, and Qm = 0.3. 

^ main subcluster only 

^ no lensing data available; estimate from the distribution of brightest galaxies 


TABLE I. List of merging galaxy clusters, with their luminosity distance di, redshift z, angular dark matter separation A@dm 
( as estimated from the distance between the two most prominent maxima in the weak lensing map), mass M 200 , and radius 
r 2 oo. The Bullet and Coma Clusters are the two benchmark clusters in this paper. M 200 and r 2 oo should be thought of as 
rough estimates for the mass and radius of the cluster, as the determination of these parameters assumes spherical symmetry 
that merging clusters do not have. We refer the reader to the references for details of the uncertainties and assumptions that 
enter into these parameters. 


We will show how the weak lensing map and the X-ray 
luminosity map, representative of the DM density and 
ICM density spatial distributions, respectively, can be 
used to improve indirect detection studies. While the ba¬ 
sic idea is simple—cluster mergers separate DM and ICM 
so that they can be distinguished—there is still a consid¬ 
erable amount of overlap between the distributions, and 
some statistical techniques are necessary to fully utilize 
the spatial information. To that end, we develop three 
major procedures in this paper: 

• Method A: If an excess of photons is observed in a 
merging cluster, the spatial distribution of the sig¬ 
nal can be used to discriminate between DM decay 
and ICM emission (or any other pair of known spa¬ 
tial distributions) as the origin of the excess. We 
develop a quantitative procedure for this discrimi¬ 
nation and show that it is quite powerful in realistic 
scenarios (e.g. the 3.55 keV line). 

• Method B: We develop a procedure to charac¬ 
terize the spatial distributions of potential DM ex¬ 
cesses by fitting to a combination of multiple spatial 
templates. This procedure is useful when the ex¬ 
cess could be from a combination of sources, and 
we specifically use it to allow for uncertainty in the 
spatial distribution of annihilating DM. This proce¬ 
dure generalizes Method A and is equally powerful. 

• Method C: Even if no definitive excess is ob¬ 
served, one can reweight the observed photons 
based on their spatial position to either enhance po¬ 
tential spectral features or to improve the exclusion 
limit (on, e.g., the DM annihilation cross-section or 
lifetime). Whereas previous studies [231 02 ] have 


masked out regions with low expected signal-to- 
background, reweighting allows all of the photons 
to contribute in an optimal way, improving the ex¬ 
pected discovery reach for DM in merging clusters. 

These procedures only depend on the spatial distribu¬ 
tions of background and signal and the angular reso¬ 
lutions of the observing instrument and lensing map, 
and they apply equally well to both broad excesses and 
monochromatic lines. 

While lensing maps do have uncertainties (a dramatic 
example is the “dark peaks” issue in A520 [35l |36j) we 
have attempted to make our methods robust against such 
uncertainties by using the entire lensing map rather than 
one or more individual peaks. We analyze the robust¬ 
ness of our procedure against lensing map uncertainties 
in Sec. |II A| and conclude that they are unimportant for 
the small signals typical of potential DM excesses. 

Prior studies have used spatial information as a diag¬ 
nostic for potential DM signals, but always in the GC 
or non-merging clusters, and based on expectations from 
simulations rather than measured distributions [451146] . 
Recently, a spatial likelihood method similar to one we 
describe in Sec. m was used in an attempt to charac¬ 
terize the origin of the 3.55 keV excess in the GC and in 
the Perseus cluster [47]. We advocate instead perform¬ 
ing these spatial tests on merging clusters (especially the 
nearby Coma Cluster), in which the spatial distribution 
of DM is both measured (by weak lensing) and quite 
different from the distribution of potential backgrounds. 
Similarly, Ref. 05 ] proposed a method to distinguish DM 
decays or annihilations into gamma rays from the cosmic 
ray-induced background in non-merging clusters. Their 
method is based on comparing the slope of the signal 
distribution with various models of DM and ICM spatial 
profiles, whereas our methods are based on the measured 
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surface mass density and use all of the available spatial 
information in an optimal way. 

Spatial information has also been used in an attempt 
to strengthen limits and boost discovery potential: the 
weak leasing map and DM-ICM separation of the Bullet 
Cluster were used to set limits on sterile neutrinos [iH] , 
However, their procedure actually resulted in a less strin¬ 
gent limit overall. Our reweighting procedure, described 
in Sec. |II C| and applied to the Bullet Cluster in Sec. |IV C] 
takes all of the DM into account and can strengthen 
the Bullet Cluster limit on the lifetime of a sterile neu¬ 
trino by a factor of ^1.2. An analogous reweighting pro¬ 
cedure, based on both spatial information and overall 
signal-to-background ratio in a stack of dSph galaxies 
was performed in Ref. [49] to search for DM annihilating 
to gamma rays. By contrast, our reweighting procedure 
(Method C), which is based on a measured weak leasing 
map, cannot directly improve searches for annihilating 
DM. However, once observed, an annihilating DM signal 
could be tested by our template fit (Method B). 

We describe our statistical methods A, B, and C in 
Sections H A[ HB and jll C[ respectively. In Section O 


we describe the range of applicability of our methods to 
indirect detection experiments in various energy ranges. 
In Sections jlV A[ |IVB[ and |IV C[ we demonstrate the 
power of our methods in a few realistic scenarios. Finally, 
we discuss the implications of the techniques proposed 
here for future indirect detection studies in Sec. B Ap¬ 
pendix]^ details an extension to our statistical reweight¬ 
ing procedure for optimizing conservative (as opposed to 
background-subtracted) limits. 


II. STATISTICAL PROCEDURE 


A. Spatial tests 


In this section, we study to what extent we can use the 
spatial information of a signal from a galaxy cluster to 
test for consistency with different spatial distributions. 
We will first consider the simple case of distinguishing 
between two dehnite spatial distributions (leaving a gen¬ 
eralization to more than two to Sec. HB[). This technique 


can be used to differentiate between two possible signal 
origins with known spatial distributions, such as DM de¬ 
cay (true signal) vs. ICM emission (background). 

Concretely, let us assume the signal sits on top of a 
large background with normalized, fractional spatial dis¬ 
tribution bi, where i labels spatial bins; e.g., a spatially 
uniform distribution with k number of bins would have 
bi = 1/k for all i. For simplicity, we take the photon 
counts in each spatial bin i to be independent random 
variables, which is a reasonable assumption as long as 
the bin size is greater than or equal to the angular reso¬ 
lution of the observing instrument. 

Suppose an excess of photons is observed in a particu¬ 
lar energy band, with a significance of s number of stan¬ 
dard deviations; assuming large statistics, this amounts 


to an observation of N+s'/N photons where only N were 
expected. If the excess of photons is due to a component 
with a fractional spatial distribution fi (with /, = !), 
then the data—given by the number of counts Xi in each 
bin—should be consistent with being a Poisson random 
variable: 


Uq: x^'^ Pois [Nbi + sVN , (I) 

a scenario which we call the null hypothesis Hq. The al¬ 
ternative hypothesis Hi is that the excess follows some 
other spatial distribution gi, in which case the data 
should be distributed according to 

Hi : Xi ^ Pois (^Nbi + sVNgi^ , (2) 


where J2i9i = 1- 

We devise a simple statistical test to distinguish be¬ 
tween these two hypotheses. As neither Ho nor Hi have 
unknown parameters, the Neyman-Pearson lemma states 
that the most powerful statistical test to distinguish be¬ 
tween them is the likelihood ratio test, which rejects Ho 
in favor of Hi when 


A(f) = F 


/Prob(af|Hi)\ 

VProb(f|Ho)y 


> A* 


(3) 


where F is any monotonically increasing function, 
and A* is determined by the test size a via 
Prob (A(a:) > A*|Ho) = a. In other words, the Neyman- 
Pearson lemma ensures that A(x) (or any monotonic 
function of the likelihood ratio) is the variable for which 
the pdfs under the two hypotheses Ho and Hi have the 
minimal amount of overlap, and therefore provides the 
strongest statistical discrimination between the two sce¬ 
narios. 

Assuming sufficiently large statistic^ in each spatially 
resolved bin —Nbi ^ I for each i —we can treat the Pois¬ 
son distributions in Eqs. Q & ([^ as Gaussians, and the 
most powerful discriminant takes the simple form 

Mi)" 2^175 W 

i * 

where we took F in Eq. ([^ to be a rescaled logarithmic 
function to simplify future expressions. The probability 
distribution of A(a:) is a weighted non-central distri¬ 
bution. There is no simple closed-form expression for this 
distribution, but each term in Eq. Q rapidly approaches 
a Gaussian for large (xi)j^and the sum over spatial bins 


^ In particular, we do not require the number of signal counts per 
bin, fi or sy/Ngi, to be large. We discuss the practicality 
of this assumption further in Sec. |III[ 

^ The skewness of each individual term in Eq. is parametrically 
suppressed as 0{Nbi)~^^‘^, and the excess kurtosis as 0{Nbi)~^. 
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i makes the pdf of A converge even faster to a normal dis¬ 
tribution (by the central limit theorernV We will treat it 
as a Gaussian in the rest of the paperjj 

Let us now extract the behavior of the A test statistic of 


Eq. 0 in the large-limit. Assuming for the moment 
that fi, and gi are perfectly known, we have for the 
expectation values and variance of A: 


(A|Hi) - (A|Ho) 


Var(A|Ho,i) 


1 9i ~ fi 

2iV3/2 2 ^ 

i—1 ^ 

1 (ffz ~ fi)"^ 

4^3 ^4 

i=l » 


(Nb, + sy/Ngif - 


Var(x2) 



4Ar3 

2=1 


{Nb, + 


2 = 1 
k 

E 


jgt - fi)'^ 

h 

{9^ - fif 

h 


T^s, 

( 5 ) 

rp‘2 

( 6 ) 


Because p(A) is well-approximated by a normal distri¬ 
bution, the expected significance of the likelihood ratio 
test can be expressed as the number of standard devia¬ 
tions: 


^ ^ (A|Hi) - (A|Ho) ^ 


v/Var(A|Ho) 


T = 




i9i - fi)" 


( 7 ) 


In the case of 0(1) separation —\gi — fi\ ^ 0{l/k )—the 
overlap factor T can be close to unity (or even exceed it). 
In addition to being the statistically most powerful dis¬ 
criminant, the test statistic A has the desirable property 
of being insensitive to the number of background pho¬ 
tons N, and the details of spatial binning, as long as the 
gross substructure of fi and gt is resolved. We shall see in 


Section IV A that A is a powerful discriminant between 
Hq and "Hi in realistic scenarios such as discriminating 
between DM and ICM origins of an X-ray excess. 

So far, we have only considered statistical contribu¬ 
tions to the variance of Xi (and thus A), ignoring possible 
systematic uncertainty in the spatial distributions fi, 
and gi- To leading order in the number of photons N, the 
statistical variance of the counts in each bin is the same 
for both Ho and Hi: Varstat(a:i) — Nbi. The system¬ 
atic contribution to the variance of Xi from an imperfect 
knowledge of bi and fi given the null hypothesis Ho is: 


Varsys(xi) ~ N^Yar{bi) -t s^lVVar(/j), 


( 8 ) 


with a similar expression (with fi ^ gi) for the alterna¬ 
tive Hi- In principle, this uncertainty should be incor¬ 
porated in the formula of Eq. Q: a larger variance in 
the pdf of A given Ho would degrade the expected power 
of the test (to exclude "Ho), while a larger variance of A 


^ We have checked via Monte Carlo that the formula in Eq. 0 
for the statistical significance s is a good approximation up to 
very large values of s for all cases of interest. In the example 
considered in Sec. uni the fractional error on s is < 1% for 
s < 8. 


given Hi implies a larger spread in the expected signifi¬ 
cance s but with the same expected value (s). However, 
the first term in Eq. Q is smaller than the statistical 
variance Nbi given a satisfactory background model bi 
(which, for example, can be measured from the side en¬ 
ergy bands). The second term in Eq. 0 can dominate 
over the statistical variance, but only for very large sig¬ 
nals s > ^/biJYax{Ji) ^ SNR(/i)-\/fc, and similarly for 
gi- For example, the leasing maps in Fig. have an aver¬ 
age per-bin signal-to-noise ratio (SNR(Ki)) greater than 
2, for fc = 50 number of bins. The uncertainty in the 
leasing map would thus only become important for sig¬ 
nals with significance greater than about s ~ 15 in this 
case. 


B. Spatial template fits 

The hypothesis test procedure described above is opti¬ 
mal (by the Neyman-Pearson lemma) for discriminating 
between two hypotheses with no free parameters. In a 
more realistic scenario there may be more than two spa¬ 
tial distributions of interest, e.g. the temperature, num¬ 
ber density of galaxies, etc., in addition to the mass and 
ICM distributions considered above. More crucially for 
our discussion, the spatial prohle of annihilating DM 
in galaxy clusters is not well known, but is likely de¬ 
scribed by a mixture of a sharply peaked core compo¬ 
nent that follows the mass density distribution raised to 
some power, and a substructure component with a more 
uniform distribution, as we will discuss further in Sec¬ 
tions |III| and |IVB[ In these circumstances, one way to 
proceed is by fitting to a linear interpolation of the rel¬ 
evant spatial templates. As above, we will work under 
the assumption that the errors on all bins can be consid¬ 
ered to be Gaussian; in that case the maximum-likelihood 
estimate of the interpolation parameters coincides with 
the minimum-y2 estimate. Additionally, the covariance 
matrix of the best-fit parameters (and, from that, the 
confidence intervals) can be extracted from the second 
derivatives of the minimum. 

Although the Neyman-Pearson lemma does not apply 
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to this procedure, we will show that in particular cases 
the fitting procedure reproduces the results of the hy¬ 
pothesis test exactly. As a result, we expect that the 
fitting procedure adequately captures all of the spatial 
information available in the signal. 

First, we describe the general case with an arbitrary 
number of candidate spatial templates, generalizing the 
two distributions / and g from above. Suppose that we 
consider fitting an s sigma excess (atop N background 
photons) to n spatial templates /“, a = 1.. .n. (We al¬ 
ways take spatial templates to be normalized to unit sum, 
fi = !■) Then we have n parameters 6a representing 
the fractional contribution of map /“ to the excess. In 
other words, our model function for the distribution of 
counts Xi (where i labels the spatial bin) is 


g^{{0a}) = (x.|{0a}) = Nb, + sVN{6afn- (9) 

We will constrain = 1, so that only n — 1 of the 

parameters are free. 

The best fit parameters da are those that minimize the 
chi-squared: 


X^{{0a}) 


(xi — gii{6a}))^ 


( 10 ) 


and the (n — 1 ) x (n — 1 ) covariance matrix Vab = 
Cov{0a, 0b), with a, 6 = 1 ,..., n— 1 , is calculated through 
its inverse as: 


(Vab)-^ = 


IdWiOa}) 


2 d0ad0b 


9^=e 




(/r - /r)(/f - m 


( 11 ) 

( 12 ) 


where the second line follows from taking Nbi large. 

The la confidence interval for 0a is 0a ± Y^Var(0a); 
Var(do) = Vaa for a = l...n — 1, and Var(0„) = 
Sa fcii lafc- Note that Var( 0 a) depends on all of the /“, 
a = 1... n, through the matrix inverse. The covariance 
matrix is still determined by geometric quantities only: 
s^Vab depends only on the spatial distributions, not their 
relative amplitudes. However, the degree to which a par¬ 
ticular hypothesis (characterized by 0 °) is excluded does 
depend on the best-fit amplitudes according to the for¬ 
mula 


S=^{dl-0a)V-,\0l-0b) 


(13) 


where s is the number of standard deviations to which 
the hypothesis is excluded. Taken as a function of 
0a, Eq. (13) defines the s-sigma error ellipsoid of the fit 


parameters. Since (x s^, s (x s exactly as in Eq. ([(^. 

Now we consider the strength of the fitting procedure 
applied to two spatial distributions and show that we 
reproduce the result of the optimal Neyman-Pearson hy¬ 
pothesis test. Eor the case of two spatial distributions, 
the fit procedure provides a slight generalization to the 


hypothesis testing procedure by allowing for a single free 
parameter 0 describing the fraction of the excess origi¬ 
nating from the spatial distribution = 5 , while 1 — d is 
the fraction of the excess originating from f^ = f. In this 
case, Eq. (12) reduces from a matrix to a single number: 


(Var(0))-i~s2^E^^E. 


(14) 


The expected best-ht parameter is 0 = 0 for the gas-only 
scenario (^. Eq. §) and 0=1 for the DM-only scenario 
(cf. Eq. (0). In either case, the expected significance 
of the result (exclusion of the alternative hypothesis) is 
simply 


s ~ 






and we recover the result of Eq. ([7]): s = Ts. 

In fact, the fitting procedure is optimal (i.e. equivalent 
to the Neyman-Pearson test) in a generalized case where 
we compare one distribution (WLOG /"■) to an arbitrary 
mixture of the rest. With null hypothesis 0° = 1 and 
expected best-fit point 0„ = 0, Eq. (13) becomes 


\ 


E 


(^a/f - /r )2 


(16) 


Since 6a ff is a properly normalized spatial template, we 
recognize Eq. ([^ with the best-fit signal distribution in 
place of gi] in other words, the fitting procedure repro¬ 
duces the optimal result of the hypothesis test despite 
not knowing the correct signal distribution in advance. 


C. Spatial reweighting 

In the absence of an obvious excess in the spectrum 
of photons from a galaxy cluster, spatial information of 
the DM distribution and the background distribution can 
still be useful to enhance a potential signal relative to 
the background by spatially reweighting the photons. A 
reweighting procedure that boosts any signal correlated 
with the dark matter distribution would aid in the iden- 
tihcation of potential features and strengthen DM ex¬ 
clusion limits if none are found. In Sec. |IV C[ we will 
demonstrate the approach outlined below on a realistic 
(though hypothetical) example of an indirect detection 
analysis. 

Suppose a total number of photons B = X)^=i 
are observed from a cluster, with Bi counts in each 
spatial bin i, and that a DM signal would amount to 
S = photons, with Si ex gi (the fractional spa¬ 

tial distribution of the signal). If the total photon count 
is dominated by the background in each bin, we expect 
Bi oc bi (the fractional spatial distribution of the back¬ 
ground). Our goal is to determine the optimal weights 
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Wi such that S = WiSi and B = WiBi yield the 
maximum statistical significance. If there is a component 
of the data which follows the DM signal distribution, it 
will be enhanced relative to the background with appro¬ 
priately chosen weights. 

To maximize the statistical significance of a signal, we 
choose weights Wi to maximize 


S({wi}) 


{s) ^ Ez wjS, 

VVarB \/E» wfBi ’ 


(17) 


the expected number of standard deviations of the excess 
above the background. Eq. 0 applies under the as¬ 
sumptions that the dominant source of error is statistical 
and that each Bi is large enough to be an approximately 
Gaussian random variable. The target function I]({wj}) 
is invariant under rescaling the weights, so we choose the 
normalization condition E^i Wi = k so that the average 
of the weights Wi is unity. Maximizing £({ 11 ;^}) subject 
to this constraint gives the optimal weights: 



(18) 


with c = fc/(E?=i Si/Bi) chosen to satisfy the normal¬ 
ization condition. Crucially, the w* are insensitive to 
uniformly rescaling Si (as any such rescaling will be ab¬ 
sorbed into c), so the total flux produced by DM need 
not be known a priori to determine the optimal weights. 

The boost in significance provided by this reweighting, 
relative to a simple integration of all bins {wi = 1 for all 
j), is given by 


R = 


s({<})_ 




E&- (19) 


i=l 


The enhancement factor R is insensitive to rescalings of 
both Si and Bi] like T (the strength of the hypothesis test 
procedure in Sec. IIA), i? is a constant that depends only 
on the spatial distributions of signal and background. 
In Sec. IV C[ we discuss the practical application of this 
technique, including typical values of R for several sce¬ 
narios. 

Spatial reweighting can also strengthen the exclusion 
limits derived from data in the absence of a signal. For 
limits derived using a background-subtraction procedure, 
the optimal weights are also given by Eq. (18). For con¬ 


servative limits that do not perform background subtrac¬ 
tion, a different weighting procedure is optimal; see Ap¬ 
pendix 1 ^ 


III. APPLICATION TO DARK MATTER 
SEARCHES 

Having developed the statistical underpinnings of our 
methods in Sec. |Tlj we now discuss their applicability to 
indirect dark matter searches in galaxy clusters. Three 


criteria need to be simultaneously satisfied for our meth¬ 
ods to be relevant: ( 1 ) applicability in well-motivated 
photon energy ranges, ( 2 ) detectable and resolvable flux 
from galaxy clusters, and (3) sufficient knowledge of the 
spatial distributions of the dominant background and 
possible dark matter signals. In this section, we show 
that these criteria are indeed fulfilled. 

The photons produced by DM decay or annihilations 
typically have energy on the same order as DM parti¬ 
cle mass, which for fermions must be above a keV due 
to phase space |50j and structure formation m con¬ 
straints. Therefore, the X-ray and gamma-ray regimes 
are the most interesting energy ranges to look for pho¬ 
tons from fermionic DM. Bosonic DM may be much 
lighter than a keV [S2] but indirect detection in the op¬ 
tical, infrared, and radio bands is subject to large astro- 
physical backgrounds; the X-ray and gamma-ray back¬ 
grounds are somewhat cleaner and their spatial distribu¬ 
tion is understood. Moreover, a DM particle (bosonic or 
fermionic) with a mass near the electroweak scale and a 
self-annihilation cross section of {av) ^ 3 • 10“^® cm^ s“^ 
is thermally produced in the early universe with the cor¬ 
rect DM relic abundance—the so-called WIMP miracle— 
motivating searches for annihilations into gamma rays. 

The second criterion for the applicability of our meth¬ 
ods is that the DM-induced photon flux from galaxy clus¬ 
ters should be observably large, and resolvable by the 
observing instrument. While our Galactic Genter at 8 
kpc is the brightest object in the sky in terms of pho¬ 
tons produced in any DM model, dwarf spheroidals and 
galaxy clusters also are interesting targets. Clusters are 
much farther away but make up for it with their large 
overall mass (see Table |l|, large DM fraction of over¬ 
all mass [53], potentially large substructure boosts for 
annihilations (see below), and the fact that the astro- 
physical backgrounds emanating from the galaxies in the 
clusters are correspondingly lower due to their large dis¬ 
tances from Earth. 

In the X-ray regime, galaxy clusters have been used as 
potential sources of photons from decaying DM, leading 
to limits competitive with those of other astrophysical 
targets |23||51|, or, in the case of the 3.55 keV line m 
a potential signal consistent with previous upper limits. 
Soft X-ray observatories also have excellent angular reso¬ 
lution, more than sufficient to resolve the gross substruc¬ 
ture of galaxy clusters, as shown in Fig. For hard 
X-rays and soft gamma rays (10 keV ^ ^ 10 GeV), 

observations out of the galactic plane set stringent lim¬ 
its on models of decaying and annihilating dark matter 
(see |3| for a comprehensive review). However, obser¬ 
vatories in this energy regime have yet to achieve the 
required angular resolving power to discern features of 
galaxy clusters (see Fig.[^. The same lack of angular res¬ 
olution will make it difficult for e.g. Fermi-LAT to check 
the Galactic excess at 1-10 GeV first reported in [10] 
with merging cluster targets. More energetic gamma-rays 
(> 100 GeV) can be sufficiently resolved by Fermi-LAT 
and air Gherenkov telescopes, as shown in Fig. [^ 

















Cherenkov telescope observations by MAGIC in the 
Perseus Cluster [55] , VERITAS in the Coma Cluster |56| , 
and HESS in the Eornax Cluster m have set limits on 
the annihilation cross section at the (av) ~ 10“^^ — 
10“^^ cm^ s“^ level in the TeV mass range, while Eermi- 
LAT observations in nearby clusters [58H60] set limits 
down to 10“^^ cm^ s“^ in the 10 GeV range. These lim¬ 
its conservatively assume a smooth NEW profile, while 
baryonic contraction EH ED, DM substructure [HHHHH] 
and Sommerfeld enhancement may substantially boost 
the DM flux to yield sensitivity to thermal cross sections 
of 3 • 10“^® cm^ s“^. For example, baryonic contraction 
is believed to boost the annihilation flux in Fornax by 
up to a factor of 2-6 [60] , while the “dumpiness” of DM 
may boost the potential annihilation flux by up to a fac¬ 
tor of 10^ in a Coma-like cluster [55]. Future Cherenkov 
telescopes such as the Cherenkov Telescope Array (CTA) 
may be able to reach sensitivity to thermal annihilation 
cross sections in galaxy clusters with sufficient observa¬ 
tion time and a high substructure boost factor [51157]. 
Clusters are also promising targets for DM decays into 
gamma rays, as considered in Refs. [5511551155] . 

We note here that for our methods to be applied ex¬ 
actly as described in Sec. [^ the total number of counts 
per bin must be large enough that Gaussian statistics is 
a good approximation. Though the number of counts per 
bin depends on many experimental details (including en¬ 
ergy resolution and binning, exposure time, and distance 
to source), the Gaussian approximation is typically jus- 
tihed for X-ray telescopes (where signals are typically on 
top of a large continuum background from the ICM) and 
for air Gherenkov telescopes (which see a large atmo¬ 
spheric background). It may not be applicable for space- 
based gamma-ray observatories such as Fermi-LAT. In 
the case of small statistics, our methods can still be ap¬ 
plied but their power will now depend on the number 
of counts per bin, making a model-independent predic¬ 
tion of their utility problematic. A full analysis of the 
small-statistics case is beyond the scope of our paper. 

Finally, the third criterion and crux for our methods is 
prior knowledge about the spatial distribution of signal 
and background. In galaxy clusters, weak gravitational 
leasing maps provide substantial information about the 
morphology of the dark matter distribution (and thus 
the signal), while the spatial distributions of the dom¬ 
inant backgrounds are extended, mostly measurable in 
(X-ray) side energy bands, and, crucially, partially non¬ 
overlapping with the signal in merging galaxy clusters. 
This is in contrast to the Galactic Genter, which has 
large uncertainties on the spatial distribution of both sig¬ 
nal and background. As we will show in Sec. |IV[ spatial 
information provides a powerful handle on even modest 
photon excesses from clusters. 

If an excess flux of photons were coming from dark 
matter decays, their spatial distribution in terms of the 
coordinates right ascension a and declination S on the 
celestial sphere would have to be correlated with the 
integral along the line of sight of the DM density dis- 



FIG. 2. (Color online) Angular resolution as a function of 
photon energy Ej for current (solid) and planned (dashed) 
X-ray and gamma-ray observatories. The light and dark gray 
bands are between 1/2 and 1/4 of the DM separation (as 
quoted in Table [^ in the Bullet and Coma Clusters, respec¬ 
tively. Our methods have reduced power inside the band, 
but apply fully below. Instrumental resolution is quantified 
as the estimated FWHM of X-ray satellites [70], CRIPS | 71| . 
INTECRAL-IBIS [72], COMPTEL [73|, ECRET [74], and 
gamma-ray telescopes m- The Hofmann limit (dotted) is the 
theoretical lower bound on the resolution of an air Cherenkov 
telescope [75] . 

tribution: oc pouia, S, 1) dl. The right- 

hand side of this equation is very nearly the surface 
mass density k measured in weak gravitational leasing: 
K(a, S) (X Pm(q:, <5, 1) dl, where the integral is over the 
total matter contribution pM = Pdm+Pvi^ which is dom¬ 
inated by dark matter so that pu — PdmIj Weak leasing 
maps such as those depicted in Fig. [^ thus provide an 
excellent model of the expected spatial distribution of a 
DM decay photon signal. For the Bullet Cluster, we used 
weak lensing data described in Ref. [22] publicly available 
at Ref. [55]; the lensing map of the Coma Cluster is the 
one used in Ref. [55] 

The spatial distribution of photons coming from an¬ 
nihilations is proportional to the line-of-sight integral of 
the square of the dark matter distribution: $“''(«, d) oc 
jy^^p^yi{a,5,l) dl. There is no other direct measure of 
this quantity, although a reasonable proxy for this ob¬ 
servable can be constructed based on the weak lensing 
map. The DM halo density of any structure is expected 
to fall off as the cube of the distance from the center of 
the halo pdm(?') cx r~^ at relatively large radius, imply¬ 
ing that p-Dmir) dl oc and /Odm(?')^ oc 
if pdm(^) is smooth. Hence the annihilation flux from 
the smooth DM halo should be tightly correlated with 
the function $“"(a, d) oc K(a, 5)®/^, i.e. a more concen- 


^ One has to take into account a convolution with the resolution of 
both the observing instrument and the weak lensing map mea¬ 
surements to match with k. 

® We thank Nobuhiro Okabe for sharing Coma lensing data. 
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trated version of the weak leasing map. However, the 
dark matter distribution is not expected to be smooth; 
in particular, it may be more “cuspy” than and/or 
“clumpy”, meaning that the quantity {p^y ^/is 
probably much larger than unity on resolvable scales, in¬ 
creasing the DM annihilation flux. Baryonic infall may 
make the DM halo core more cuspy [50] . 


Self-bound DM subhalos are also expected to signifi¬ 
cantly boost the annihilation flux. Simulations of cold 
DM halos on the cluster scale [Ml ESI El] with resolu¬ 
tions down to ~ IO^Mq show that these resolved sub¬ 
structures contribute roughly twice the annihilation flux 
of the NFW-like core. However, substructure may exist 
down to scales of 10“® to depending on the 

nature of the DM and perturbation spectrum, in which 
case the estimate of the flux from the unresolved sub¬ 
structure ranges from ~ 50 to ~ 1000 times the core 
flux. The spatial distribution of the annihilation flux 
from substructure is expected to be approximately uni¬ 
form throughout the cluster up to some cutoff radius, 
after which it falls off as r“^, similar to the expected flux 
from DM decays [55]. In this paper, we construct tem¬ 
plates for the core and substructure distributions from 
the weak leasing map. We then allow the relative con¬ 
tribution of core and substructure to vary in accordance 
with our fitting procedure. (The details are addressed 
in Sec. jlVB ) We present our results as estimates of the 
strength of our procedure as applied to annihilating DM, 
and hope that future simulations will shed more light on 
the expected substructure in merging clusters. 


For the spatial distribution of the ICM plasma, we use 
public data from X-ray observatories: XMM-Newton for 
the Bullet Cluster [29], and ROSAT for the Coma Clus¬ 
ter |5D]. The flux of X-rays from the ICM plasma is 
dominated by thermal bremsstrahlung with only a weak 
dependence on temperature; hence the intensity of X- 
rays is a measure of the (squared) density of the plasma. 
The data amounts to a direct measurement of the back¬ 
ground in the X-ray energy range, and we also use it 
as a proxy for the distribution of other features corre¬ 
lated with the gas emission. In the gamma ray energy 
range, only the (spatially uniform) diffuse extragalactic 
background has so far been observed in clusters. (Other 
gamma-ray backgrounds, originating within the Galaxy, 
may also contribute to the observed background; for our 
purposes, we will assume that these are subdominant to 
the spatially uniform part.) However, the dominant non- 
DM source of gamma rays from clusters is expected to be 
produced by scattering of cosmic rays on the ICM. The 
cosmic rays (CR) are well-confined by the ICM’s mag¬ 
netic field with a long cooling time, resulting in a flux 
proportional to ucruicm — ’^icm therefore similar 
in spatial distribution to the X-rays [78] . 


The normalized spatial distributions that we will use in 
this paper are displayed in Fig. The ICM and DM de¬ 
cay distributions are measured data (X-ray flux and weak 
leasing convergence, respectively) as discussed above; the 
DM core and DM substructure distributions are approxi¬ 


mations derived by combining the weak leasing map and 
phenomenological models from simulations. The exact 
constructions are discussed further in Sec. IIVBI 


IV. RESULTS 

Having developed our statistical techniques in Sec. jH] 
and the expected spatial distributions of signal and back¬ 
ground in Sec. |HI[ we now move on to cmantify the power 
of these methods in realistic situations lH We will discuss 
the application of (A.) the spatial test procedure to an 
X-ray line, (B.) the template fitting procedure to anni¬ 
hilating WIMP DM in gamma rays, and (C.) boosting a 
small X-ray excess using spatial reweighting. 


A. Spatial test: Sterile neutrino decay 


We now quantify the power of the hypothesis test out¬ 
lined in Sec. |H A| Suppose an indirect detection experi¬ 
ment sees an anomalously large number of photons, and 
that this excess has a significance of s sigma. We would 
like to know whether this excess of photons has a spatial 
distribution that is consistent either with a DM origin (in 
which it should be correlated in some way with the weak 
leasing map), or with an unmodeled component of the 
ICM. Specifically, we want to establish how well one can 
distinguish an excess of s sigma significance with spatial 
distribution from one with spatial distribution on 
top of a large background with spatial distribution bi. 
(As in Sec. HA spatial bins are labeled by i and the dis¬ 


tributions are normalized so that ft, gi,bi = 1.) As a 
demonstrative example, we will use the 3.55 keV line m 
to show how the method works in practice. 

For this particular anomaly, the main background is 
composed of the X-ray emission of the ICM, of which the 
spatial distribution bi can be measured accurately in side 
energy bands. The most conservative hypothesis is that 
the observed excess is an unmodeled emission line of the 
ICM, in which case it should spatially follow the main 
background, i.e. fi ~ bi. An exciting possibility is that 
the excess is due to a dark matter decay such as that of 
a sterile neutrino into an X-ray photon and a neutrino, 
in which the signal should be strongly correlated with 
the weak leasing convergence Ki. Therefore, our alterna¬ 
tive hypothesis is that the excess has a spatial distribu¬ 
tion Qi oc Ki- In Sec. |H A[ we showed that the optimal 
way to discriminate between these two possibilities is to 
compute the quantity A = (2fV^/^)“^ ^-(g^ — /i)a:f/ 6 f, 
where Xi is the observed number of counts in each spa¬ 
tial bin. Since the expected value of A is different for 


^ Even though we illustrate our methods on real X-ray data, we 
refrain from doing a careful analysis; we leave this to the ex¬ 
perts with a better understanding of background modeling, point 
source subtraction, and instrumental response/calibration. 
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Bullet - ICM Bullet - DM Decay Bullet - DM Substructure Bullet - DM Core 



FIG. 3. (Color online) Normalized spatial distributions of signals originating from gas, decaying DM, and annihilating DM 
(substructure and core components) in the Bullet (top) and Coma (bottom) Clusters. 


the two possibilities we are considering, the value of A 
computed from the data can tell us which distribution is 
a better match. If the ICM is the source of the excess 
(e.g. an ICM emission line), we would expect (A) ~ 0; if 
the excess originates from decaying DM, then (A) ~ T^s, 
where ~ fi)^/bi is determined by the sepa¬ 

ration of the spatial distributions of the ICM and DM. 
The standard deviation of A under either hypothesis is 
T, so the two hypotheses can be distinguished at s = Ts 
sigma. 

Suppose an X-ray telescope with effective area of about 
400 cm^ had an exposure of a few time s 10® s to the 
Coma Cluster region depicted in Fig. IbP Such an expo¬ 
sure would yield a very large number of background ICM 
photons in the energy band around 3.55 keV {N > 10®, 
depending on the energy resolution of the instrument), as 
well as potentially a few hundred “excess” photons from 
the unidentified emission line reported in m- For defi¬ 
niteness, we assume that this experiment led to an s = 5 
sigma detection of the emission line at 3.55 keV in the 
energy spectrum, not unlike what was reported in 
Spatial binning into square bins of 4', the FWHM an¬ 
gular resolution of the weak lensing map, as depicted in 
Fig. lb (for a total number of spatial bins k = 50) satisfies 
the necessary criteria for our procedure to be effective. 
The number of counts Xi in each bin is sufficiently large, 
Nbi 1, that they are appropriately described by Gaus¬ 
sian random variables with Vara;, ~ Nb^. Furthermore, 
the angular resolution of a typical X-ray telescope (on 
the order of 10") is much smaller than 4', so we can treat 


® For calculations in this section, we use ROSAT data EQl which 
does not cover the 3.55 keV energy band. XMM-Newton obser¬ 
vations of Coma exist but require combining several exposures. 
Since we do not expect the spatial distribution of the X-ray back¬ 
ground to change significantly from 2-3.55 keV, we chose the 
ROSAT data for simplicity. Combining the XMM-Newton expo¬ 
sures adds a technical complication but does not affect the power 
of our test procedure. 

^ We use the 3.55 keV line only as a specific case study for il¬ 
lustrating the implementation and power of our methods, and 
acknowledge the existing tension between different studies of the 
excess, as discussed in Sec. IT] 


the counts in each spatial bin as independent. 

The power of the proposed test is depicted in Fig. 
Even with this relatively coarse binning, the ICM and 
DM distributions of the Coma Cluster are well-separated, 
with discrimination factor T « 0.68. This causes the 
probability density functions (pdfs) p(A|ICM) (blue, left) 
andp(A|DM) (red, right) to differ substantially: the ex¬ 
pected value of A under the DM hypothesis is Ts « 3.4 
standard deviations removed from the expected value of 
A given the ICM hypothesis. For the case at hand, the 
above hypothesis test on the s = 5 sigma excess could 
rule out an additional ICM line contribution in favor 
of a DM decay scenario with an expected significance 
of Ts ± 1 « 3.4 ± 1 sigma. Conversely, since the vari¬ 
ances of p{A\ICM) and p{A\DM) are very nearly equal, 
a DM decay hypothesis could be ruled out in favor of 
an ordinary ICM component at the same expected sig¬ 
nificance. If the excess was less prominent, the power 
of the hypothesis test decreases also, but always accord¬ 
ing to the relation s ^ Ts: e.g. an s = 3 sigma excess 
could be tested for an ICM contribution versus DM de¬ 
cay at s « 2.0 ± 1 sigma. Hence, even moderately-sized 
excesses can be tested for consistency with a true signal 
(DM) or background (ICM) with high significance due to 
the spatial separation of DM and gas in merging galaxy 
clusters. 


A large discrimination factor T is thus all that is 
needed to distinguish a DM decay signal from an ICM 
source based on its spatial distribution. In Table |ITj 
we list the discrimination factor T for several scenarios. 
The left column is for distinguishing a DM decay signal 
from an ICM source, as discussed above. The right col¬ 
umn lists the T factors to differentiate DM decay from 
a gamma-ray source. We took the baseline gamma ray 
background bi to be uniform, and the alternative spa¬ 
tial distribution of a potential excess to be that of 
cosmic rays scattering off the ICM (and, therefore, well- 
approximated by the X-ray distribution; see Sec. Ill I. For 
both the Bullet and Coma Clusters, T is large (> 0.5) 
if the instrumental angular resolution 69 is sufficient to 
resolve the separation of the DM and gas in the cluster. 
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FIG. 4. (Color online) Probability density functions of A if an 
s = 5 sigma excess is an unmodeled emission line in the ICM 
(blue, left) or a DM decay line (red, right). An ICM back¬ 
ground of 10® photons originating from the Coma Cluster was 
injected with an s = 5 sigma statistical excess following either 
the ICM or the cluster mass density p rofi le; the photons were 
binned into square bins of 4' as in Fig. |lb[ The discrimination 
power of this test is s = Ts ± 1 « 3.4±1; in other words, 
if the excess is truly from DM, an ICM interpretation of the 
excess is expected to be excluded at 3.4 ± 1 sigma (and vice 
versa). 



X-ray 

gamma ray 

DM decay signal 

Qi 

lensing map n 

alternative 

/« 

2 

^ICM 

ncnnicM 

background 

hi 

2 

^ICM 

uniform 

59 

T 


12" 

0.79 

0.67 

Bullet Cluster 

CO 

0 

0.74 

0.62 


2.5' 

0.12 

0.11 

Coma Cluster 

4' 

0.68 

0.59 


TABLE II. List of discrimination factors T in our benchmark 
merging clusters for various angular resolutions 59 and choices 
of background (distributed like the ICM for X-rays, and uni¬ 
form for the diffuse gamma ray background). An s-sigma 
excess of photons can be tested for consistency with decaying 
DM compared to an ICM-based source (e.g. X-ray lines or 
cosmic ray-ICM scattering producing gamma rays) at s = Ts 
sigma. The necessity of sufficient angular resolution is evi¬ 
dent for the Bullet Cluster T-values; a 2.5' resolution is barely 
sufficient to resolve the two subclusters, leading to a tiny dis¬ 
crimination factor. On the other hand, once the structure is 
resolved, improving the angular resolution does not increase 
T significantly. 


B. Spatial template fit: WIMP annihilation 


To demonstrate our methods in the case of annihilating 
DM, we turn to the fitting procedure for multiple spatial 
distributions described in Sec. |IIB While the spatial dis¬ 
tribution of the decay products of DM can be measured 
directly by weak lensing, the distribution of photons pro¬ 
duced in annihilations depends on the integral dl 


which can have a drastically different shape than the sur¬ 
face mass density. For a smooth NFW profile, the surface 
mass density k falls off like r~^ at large r, while the an¬ 
nihilation signal falls off like r“®. However, substructure 
at scales smaller than what weak lensing can resolve is 
expected to contribute substantially to the annihilation 
signal; in fact, simulations indicate that it may dominate 
the total flux, as we discussed in Sec. m In the absence 
of a definitive prediction for the distribution of the an¬ 
nihilation signal in a merging cluster, we will treat the 
annihilation signal as a mixture of an NFW-like “core” 
and a flatter “substructure” component, both of which 
we will estimate from the measured surface mass density. 
Based on the large-r behavior for the NFW profile, we 
trade r for and estimate the core component as 

fcore ^ This distribution has the qualitatively cor¬ 
rect feature that DM peaks become much sharper. For 
the substructure, [5S] find that the distribution of the an¬ 
nihilation signal from substructure can be described by 
(1-b (4r/r2oo)^)~^5 i-e. roughly uniform inside r2oo/4 and 
falling off like outside. Since the distribution in terms 
of r is not very useful for a non-spherically-symmetric 
merging cluster, we again trade r for and model 

the substructure as oc (1 -I- 16 k20o/«^)~^ where K 200 
corresponds to the surface mass density at r 2 oo- Since our 
chosen region of the Coma Cluster is comparable in size 
to r 2 oo ~ 1 Mpc, we pick K 200 to be the smallest value 
along the boundary of the region. We stress that unlike 
the spatial profile of DM decay (which is directly mea¬ 
sured by gravitational lensing), there is a large systematic 
uncertainty in the spatial profile of a DM annihilation sig¬ 
nal for both the core and substructure components. The 
templates we construct from the weak lensing map are 
simply guesses that are consistent with the long-range 
behavior of DM profiles. However, we think they can 
be regarded as useful proxies: any linear combination of 
DM core and substructure components should be qual¬ 
itatively very different from the ICM flux in a merging 
cluster. Lacking any observational evidence or simula¬ 
tions of merging clusters to suggest a better approach, 
we will continue using these constructed templates with 
the caveat that our results for annihilating DM will be 
only estimates. 

We will demonstrate the discriminating power of the 
fitting procedure in the case of DM annihilating to 
gamma rays in the Coma Cluster. The basic idea is to 
find the best-fit combination of spatial templates to re¬ 
produce the excess and determine the extent to which 
various hypotheses (i.e., particular combinations of tem¬ 
plates) are excluded. Using the formalism introduced in 
Sec. |HB[ we consider htting a potential excess to the 
distributions and ^ jX-ray^ 

two free parameters, 6*core and O^uh, which are (respec¬ 
tively) the fractional contribution of the DM core and 
DM substructure to the excess. The ICM contribution is 
constrained to be 6*icm = 1 — ^core — 6*sub- We take the 
background distribution bi = uniform for gamma rays 
and the cosmic ray induced distribution to be the same as 
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the X-ray distribution, as discussed in Sec. Ill The abil¬ 
ity of the fit to distinguish between the three scenarios 
is encoded by the inverse covariance matrix of 9i — 0 core 
and 6*2 = ^sub, calculated via Eq. (Il^: 




(/f - /?)(/? - IT) 


= s 


0.811 0.259 
0.259 0.189 


( 20 ) 


From this we can determine the Icr error bars on each 
9a independent of the best-fit values: tJcore — 1.5s“^, 
tJsub — 3.1s“^, and cticm — 2.4s“^. However, the exclu¬ 
sion of the ICM-only hypothesis depends on the best-fit 
parameters via Eq. (13); when the substructure boost 
is large, 0sub —>■ 1 and s 0.43s. The inverse covari¬ 
ance matrix also defines the error ellipse in the 0 core“^sub 
plane. 

In Fig. we depict potential outcomes of our fit¬ 
ting procedure for an s = 5 sigma excess in the Coma 
Cluster. A spatial fit to an observed excess would yield 
a best-fit value anywhere in the 0 core“^sub plane with 
an error ellipse determined by the covariance matrix of 


Eq. (12). Note that the size and shape of the error el¬ 


lipse are independent of the best-fit point. The red ovals 
represent 1 - and 2 -sigma contours for a best-fit value of 
6*core = 1 and 9suh = 6*icm = 0 , i.e. consistent with an¬ 
nihilation produced mainly by an NFW-like core distri¬ 
bution. In this case, an explanation of the excess via 
cosmic rays scattering off the ICM would be excluded 
at 4.5 sigma. The blue ovals are for a best-fit value 
of dsub = 1 and 6 *core = ^ICM = 0, Consistent with a 
substructure-dominated signal, which would rule out a 
cosmic ray explanation at 2.2 sigma. 


C. Spatial reweighting: dark matter decay 



- 0.5 0.0 0.5 1.0 1.5 


^core 

FIG. 5. (Color online) The expected 1- and 2-sigma confi¬ 
dence regions for fitting the spatial profile of a 5-sigma excess 
in the Coma Cluster. 0core and 0sub are the fractions of the 
excess originating from DM annihilation in the NFW-like core 
and clumpy substructure, respectively; the remainder of the 
excess is assumed to originate from cosmic ray-ICM scatter¬ 
ing. The red ellipses (bottom right) correspond to the errors 
on the fit parameters for a perfect fit to a core-shaped DM an¬ 
nihilation flux distribution, while the blue ellipses (top left) 
are for a perfect fit to a flux dominated DM substructure. 
The triangle corresponds to the region of hypotheses that we 
are testing (although statistical fluctuations may cause the 
best-fit point to lie outside the triangle), and the green side 
(hypotenuse) is the set of DM-only hypotheses. In these co¬ 
ordinates, the shape of the confidence region (i.e. size of error 
bars) is independent of the best-fit point, but the exclusion of 
the ICM-only hypothesis does depend on the best-fit point, 
even if it lies along the DM-only line. 


We now wish to apply the spatial reweighting tech¬ 
nique described in Sec. |IIC| to physical situations. Meth¬ 
ods A and B in the previous sections were designed to test 
for a DM origin of an already-existing statistical excess. 
This procedure is optimized for “DM bump hunting”, 
i.e. to amplify putative DM-induced photon excesses in 
the spectrum. The reweighting procedure also consti¬ 
tutes a test for DM origin of an excess, though a less 
powerful one than the one we devised in Sec. |IV A[ We 
shall mostly restrict ourselves here to the case of potential 
photon excesses from DM decay, as the spatial distribu¬ 
tion from an annihilation signal is less certain. 

As detailed in Sec. mg we can reweight each pho¬ 
ton with a factor dependent on its spatial position in 
the sky. The weight factors are chosen to amplify pho¬ 
tons from DM-rich regions of a cluster and diminish those 
from background-rich regions. As long as the total num¬ 
ber of observed counts Bi in each spatial bin i is large 
enough to be approximately Gaussian, and as long as the 
dominant error is statistical, we showed that the optimal 


weights are 


w*i = c-^, ( 21 ) 

where Si is proportional to the expected number of signal 
photons and c is an arbitrary normalization factor. In 
this work, we choose c = k/{J2i^i Si/Bi) so that the 
typical value of a weight is 1 . 

We also showed that the expected boost in statistical 
significance of a true DM signal is given by 




( 22 ) 


where s (s) is the significance of the signal before (after) 
reweighting. Because R (like w^i) is insensitive to rescal¬ 
ings of Si and Bi, it is a purely geometric property of a 
particular cluster, depending only on the spatial distribu¬ 
tions of DM and background. However, to the extent that 



















13 



X-ray 

gamma ray 


0.7 

1.2 

0.8 

0.4 

0.3 

0.4 

0.8 

1.5 

1.3 

L 3 

DM decay signal Si 
background Bi 

lensiuB mao k 












2 

'^ICM 

uniform 


0.9 

0.8 

0.5 

0.3 

0.3 

0.3 

0.3 

1.5 

2.8 

1.4 



R 


0.9 

0.7 

0.6 

0.3 

0.3 

0.3 

0.4 

0.8 

2.8 

1.3 

Bullet Cluster 

1.243 

1.06 












Coma Cluster 

1.21 

1.16 

7j 

Q 

1.1 

0.9 

0.6 

0.4 

0.4 

0.4 

0.6 

1.6 

2.5 

0,9 

average value; R{E) is 

weakly energy dependent, see text 


1.3 

1.1 

0.9 

0.6 

0.4 

0.5 

1.2 

1.2 

2.1 

1.7 


boost R for DM decaying to either X-rays or gamma rays, 
in onr two representative clusters, calculated according to 
Eq. (221. A true decaying DM signal with unweighted sig¬ 


nificance s will be boosted by reweighting to a significance 
s = Rs. The regions and bins used are those of Fig. The 


bin size is 30 
Cluster. 


for the Bullet Cluster and 4' for the Coma 


the shape of the background energy spectrum varies with 
position, e.g. with temperature variations of the ICM in 
the X-ray energy band, R will vary with energy. For DM 
decaying to X-rays, the energy dependence is quite weak, 
varying from i?(l keV) « 1.3 to i?(5 keV) « 1.2 in the 
Bullet Cluster. For DM decaying to gamma rays, the 
background is uniform and independent of energy. We 
report values of R for both the Bullet and Coma Clus¬ 
ters in Table imi 

The reweighting technique is not effective for DM an¬ 
nihilating to gamma rays if the flux from the substruc¬ 
ture dominates the flux from the smooth core, since in 
that case the spatial distributions of the signal and back¬ 
ground are both very nearly uniform within the cluster 
and the spatial distribution provides no extra informa¬ 
tion. We find R < 1.05 in both the Bullet and Coma 
Clusters. The hypothesis test (A) and fitting procedure 
(B) are still effective in this scenario because the excess 
is already in hand, and we are only comparing its shape 
to an alternative (e.g. cosmic ray-ICM collisions) which 
is non-uniform. Here, the alternative is a fluctuation 
of the background (with uniform spatial distribution), 
which is practically indistinguishable from DM annihila¬ 
tion dominated by substructure flux (also nearly spatially 
uniform). 

Spatial reweighting can enhance DM spectral features 
before they are prominent enough to be considered for 
the hypothesis test of Sec. |IV A| or the fitting procedure 
of Sec. |IVB[ As a concrete example, we inject a 2.5 ct 
excess, spatially distributed like decaying DM, into the 
X-ray spectrum of the Bullet Cluster at A = 4.25 keV. 
We use the 30" spatial bins shown in Fig. la We then 


calculate the optimal weight factors according to Eq. (21) 


with Bi given by the total number of counts in a 1 keV- 
wide bin. (Although R is only weakly energy dependent, 
the weight factors themselves vary more strongly; the op¬ 
timal weights calculated at low energies are not effective 
for excesses at higher energies. However, the wide energy 
binning we used is sufficient to enhance the excess.) The 
weight factors are displayed in Fig. The unweighted 


Right Ascension a 

FIG. 6. Distribution of decaying DM signal photons in the 
Bullet Cluster (in grayscale) in 30" bins overlaid with the 
weight factors Wi calculated according to Eq. |21| The region 
and binning are identical to those of Fig. |la| The background 
is taken to be the X-ray photons in the range 3.75-4.75 keV, 
i.e. the weight factors used at the energy of our injected signal 
(cf. Fig. 0. Note that the bins with the largest weights are 
not those with the highest DM content (used in Ref. [43]), 
but those on the outskirts of the cluster where the DM has 
moved beyond the ICM and is correspondingly more pure. 


and reweighted spectra are shown in Fig.0 the injected 
DM excess at 4.25 keV, which was only 2.5cr above back¬ 
ground in the unweighted case, is visibly enhanced to a 
3(7 excess, in accordance with our calculation of i? « 1.2. 
The blue line is a simplistic estimate of the background 
contribution, calculated with by averaging side energy 
bands, meant only to aid in estimating the significance 
of the injected signal. Other l-2cr excesses visible in the 
unweighted spectrum, e.g. at 3.75 keV, are in fact sup¬ 
pressed by the reweighting. 

Spatial reweighting can also be used to strengthen 
exclusion limits by an amount proportional to R. A 
background-subtracted fc-sigma exclusion limit on some 
DM parameter a, where the received flux in the detec¬ 
tor is <i)DM oc a”, is set by enforcing that < 1 >dm produce 
no more than a k sigma excess above the observed back¬ 
ground. Since the significance of such a signal would 
be boosted by R after reweighting, the limit on a will 
be multiplied by a factor of For example, the 

flux from decaying sterile neutrino DM is linearly propor¬ 
tional to the parameter sin^ 20, and so reweighting can 
make those limits stronger by a factor of R. In Ref. |43j . 
the regions selected for their method correspond to a non- 
optimal choice of weight factors in our method. In fact, 
the i?-value calculated using these regions is less than 
1; selecting only the DM peaks gives a weaker limit on 
sin^ 29 than cutting out no region at all. This can be un¬ 
derstood from Fig. 0by noting that the weights in DM- 
rich regions are often smaller than one, while the regions 
on the outskirts of the cluster (where there is some DM 
but almost no gas) are heavily weighted. By reweighting 
the photons optimally, the limit could be improved by a 
factor of i? « 1.2 instead. 

A similar method (with different weight factors) is ap¬ 
plicable for strengthening conservative, non-background- 
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FIG. 7. (Color online) A sample energy spectrum of the Bullet Cluster with 30” binning, before (left) and after (right) 
reweighting the photons based on their spatial distribution. The blue solid line is the estimated background contribution, 
calculated from the sidebands (average of two bins on each side). A fictitious 2.5cr line excess with the spatial distribution of 
decaying DM has been injected a.t E = 4.25 keV. The i?-factor at this energy is roughly 1.2, and the significance increases 
accordingly to 2.5Ra « 3a after reweighting. The weights are calculated using 1 keV-wide bins to reduce their variance. 


subtracted limits; see Appendix [X} 


V. DISCUSSION & FUTURE OUTLOOK 


We have developed three data analysis methods for in¬ 
direct detection of dark matter in galaxy clusters, where 
the extended dark matter spatial distribution can be 
measured via weak gravitational leasing. All three meth¬ 
ods A, B, and C rely on this spatial information to en¬ 
hance the discrimination power between dark matter sig¬ 
nals and astrophysical backgrounds. 


Method A is a spatial test on a putative statistical 
excess, and assesses the compatibility of the excess 
with two known morphologies. It is especially pow¬ 
erful to determine whether an excess is consistent 


with a dark matter decay, as shown in Sec. IV A 


• Method B is a spatial fit of any number of spa¬ 
tial templates to a putative statistical excess. For 
two templates, it is equivalent to Method A, but it 
can be extended to three or more templates. It is 
most useful for distinguishing dark matter annihi¬ 
lation signals (which is expected to have a smooth 
halo component, and an extended subhalo compo¬ 
nent in galaxy clusters, with unknown ratio) from 
background, as explained in Sec. |IVB[ 


• Method C is a spatial reweighting of photons, de¬ 
signed to enlarge any statistical excess with known 
spatial distribution. It applies to both decays and 
annihilations, although it is more promising for the 


former, as discussed in Sec. IV C 


Our methods have several important advantages: 

• Robustness: In contrast to methods used in the 
Galactic Center and dwarf spheroidals, each of our 
methods are based on spatial distributions for sig¬ 
nal and background which are both measured and 
separated. For decaying DM, the spatial distribu¬ 
tion is directly measured by weak lensing; for anni¬ 
hilations, we combine the weak lensing map with 
general results from simulations. Method B ac¬ 
counts for some of the uncertainty in the annihi¬ 
lation template. The background distributions are 
also well-known and substantially different from the 
signal distributions due to the cluster merger. 

• Universality: Our results are insensitive to details 
of spatial binning (as long as the gross substructure 
of the merging cluster is resolved), energy range, 
exposure time, etc., and only depend on the intrin¬ 
sic spatial distributions of signal and background. 

• Optimality: For the questions they try to answer. 
Methods A and C are provably optimal, and use 
all of the spatial information available. Method B 
generalizes Method A and also appears to use all 
of the spatial information. 

• Wide range of applicability: Our methods ap¬ 
ply to any type of emission for which the spatial 
distributions of signal and background are known 
or can be reliably estimated. We showed how weak 
lensing maps of galaxy clusters can be used to test 
potential decaying or annihilating dark matter sig¬ 
nals in the well-motivated X-ray and gamma-ray 
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energy bands. This type of analysis can be done 
offline and with existing data. The statistical meth¬ 
ods presented here may also And a use outside the 
realm of dark matter detection, in areas where a 
weak signal is to be extracted and tested for con¬ 
sistency with an independently measured property 
of the expected signal. 

We have also showed that these methods are promising 
in realistic scenarios, including current indirect detection 
anomalies. Method A can be applied to the 3.55 keV line 
in already existing observations of the Coma Cluster. A 
longer observation of Coma leading to a 5a observation 
of the line could distinguish a decaying sterile neutrino 
from ICM emission at more than 3cr. Fermi-LAT, HESS, 
and (in the future) CTA have the angular resolution to 
test future excesses in high-energy gamma rays in the 
Coma Cluster. 

The main downside of our method is that galaxy clus¬ 
ters are fainter targets in terms of DM-induced photon 
flux than the Galactic Center (and perhaps also dwarf 
spheroidals). As a target for discovery of an anomalous 
excess, the latter is thus more interesting; to elucidate 
the DM nature of an excess, we have argued that merg¬ 
ing galaxy clusters are the most compelling targets. In 
this sense, indirect searches of DM in clusters are com¬ 
plementary to those in the Galactic Center and dwarf 
spheroidals. 

Our methods may become even better as more merging 
galaxy clusters are discovered. The well-known and very 
nearby Coma Cluster was only recently lensed and con¬ 
firmed to have undergone a signiflcant merger leading to 
DM-ICM separation, and galaxy clusters with major and 
minor mergers continue to be found |Z9|. Although Coma 
is likely the closest cluster with a substantial merging 
event, and gravitational leasing becomes prohibitively 
difficult for closer targets, future leasing studies may dis¬ 
cover interesting DM morphologies in other nearby clus¬ 
ters. New instruments with good angular resolution may 
close the hard X-ray/soft gamma ray gap (see Fig.|^ and 
allow observations of merging clusters in interesting new 
energy ranges. 

The data analysis methods outlined in this paper make 
indirect detection in galaxy clusters more attractive than 
before. In the near term, they can be applied to existing 
anomalies such as the 3.55 keV line, and provide extra 
motivation for constructing weak gravitational leasing 
maps of nearby galaxy clusters. In the future, they mo¬ 
tivate dedicating substantial instrumental exposure time 
to galaxy clusters, as well as pushing the boundaries on 
angular resolution and effective area of future X-ray and 
gamma-ray observatories. 
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Appendix A: Conservative limit reweighting 

If no excess is seen, a limit on the DM flux can be set. 
When the background is well-modeled, a background- 
subtracted limit may be a reasonable approach; spatial 
reweighting can improve such a limit only by a factor of 
R as described in Sec. |H C| and quantihed in Sec. |IV C| 
However, if the background is not well-understood, the 
most conservative way to set a limit is to constrain the 
DM flux to be less than the total observed flux (plus 
uncertainty): 


‘I’dM < $obs + (^1) 


where s is the number of sigma corresponding to the con¬ 
fidence level of the limit (typically s = 2, for a 95% con¬ 
fidence limit). 

With spatial information, we can reweight the photons 
to set a stronger limit on the DM hypothesis. In terms 
of the weighted variables S and B, the limit is set by 

K{a){S) < {B) + sVvarB where K gives the observed 
counts per unit Si in terms of the DM parameter a. For 
example, in the case of DM decay, a is some coupling in 
the model and Kip) gives the number of observed counts 
per unit surface mass density on the sky, taking Si = ki. 

Since the strongest limit is set by minimizing the lim¬ 
iting value of K{a), we choose to maximize the target 
function 


Scl({w^*}) 


(^) 

(R) -ksVVarR 
_ Yji _ 

Y.iW^B^ + s^JY~wfBi 


(A2) 


Again, the target function is unchanged by rescaling Wi 
so we impose There is no closed form 

solution for the optimal weight but it satisfies the 
equation 


We thank Masha Baryakhtar, Marusa Bradac, Kiel 
Howe, Xinlu Huang, David E. Kaplan, Jeremy Mardon, 
Ondrej Urban, and Yue Zhao for discussions. We are 


W^j = c 


( W*zBi + Sy/X), W%B, Sj _ \ 

Bj j 


(A3) 
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which can be efficiently solved by iteration if the weights 
are restricted to be nonnegative. These weights are in¬ 
sensitive to rescalings of Si but not to those of Bi] as the 
total number of counts increases, the optimal weighting 
for setting a conservative limit changes, eventually ap¬ 


proaching a Kronecker S at the location of greatest SijBi. 

The boost in limit-setting power is given by i?cL = 
ScL({w*i})/ScL({l})- Typically i?cL ^ ©(few). For 
the same exposure and binning of the Bullet Cluster dis¬ 
cussed in Section IV C[ Rcl « 3 for s = 2 and grows 
slowly with exposure time. 
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